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ABSTRACT 


It  is  very  well  known  that  Newton's  interpolation  series 

f(x)  -  H Xq )  +  (x-x0)f(xQ,x1 )  +  (x-x0)(x-x1  )f(xQ,x1,x2)  +  ••• 

simplifies  considerably  in  the  case  that  the  points  xn  *  a  +  nh  form  an 

arithmetic  progression.  Indeed,  in  this  case 

f(a,a-*h,»  ••  ,a+nft)  -  — - —  A^f(a)  . 

nlh 

It  seems  much  less  known  that  a  similar  simplification  occurs  in  the  case  when 


the  points  of  interpolation  form  a 


trie  progression.  This  paper  deals 


with  this  interpolation  problem  and  its  main  contribution  is  to  call  attention 
to  the  references  [6],  [5],  [3]  to  the  work  of  Stirling  (1730),  Schellbach 
(1864),  and  Runge  (1891),  which  seems  now  practically  forgotten.  This  work  is 
here  described  and  also  its  close  connection  with  the  elegant  algorithm  of 
Romberg  (See  [1]).  We  illustrate  these  connections  with  numerical  examples. 
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SIGNIFICANCE  AND  EXPLANATION 


"it  is  very  well  known  that  Newton's  interpolation  series  with  divided 
differences  simplifies  considerably  in  the  case  that  we  interpolate  in  the 
points  of  an  arithmetic  progression.  It  seems  much  less  known  that  a  similar 
simplification  occurs  in  the  case  when  the  points  of  interpolation  form  a 
geometric  progression.  We  describe  here  the  practically  forgotten  work  of 
Stirling  (1730),  Schellbach  (1864),  and  Range  (1891),  and  its  connection  with 
the  elegant  and  more  recent  algorithm  of  Romberg  (1955). li 
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ON  POLYNOMIAL  INTERPOLATION  IN  THE  POINTS  OF  A  GEOMETRIC  PROGRESSION , 
STIRLING,  SCHELLBACH,  RUNGE  AND  ROMBERG 

I.  J.  Schoenberg 

1*  Introduction .  It  is  very  well  known  that  Newton's  interpolation 
series 

(1.1)  f ( x)  -  f(Xg)  +  (x-Xq  )f(x0,x1 )  +  (x-x0)(x-x1)£(x0,x1,]o{)  +  ••• 
simplifies  considerably  in  the  case  that  the  points  xn  -  a  +  nh  form  an 
arithmetic  progression.  Indeed,  in  this  case  we  have 

(1.2)  f(a,a+h,««*,a+nh)  «  — - —  a"  f(a)  . 

n)hn  h 

It  seems  much  less  known  that  a  similar  simplification  occurs  in  the  case  when 
the  points  of  interpolation  form  a  geometric  progression,  at  least  this  is  not 
mentioned  in  the  standard  treatises  on  this  subject. 

The  main  contribution  of  the  present  paper  is  to  call  attention  to  the 
references  [6] ,  [5] ,  [3] ,  to  the  work  of  Stirling,  Schellbach  and  Runge,  which 
seem  now  to  be  practically  forgotten.  These  were  known  to  the  author  since 
1943.  Stirling  expands  the  function  F(z)  in  the  form 

F(z)  -  Sq  +  a.,rz  +  a2r2z  +  •••  (r  >  1)  , 

determining  the  coefficients  Sq^  ,a2,  **•  by  the  interpolation  at  z  - 
0,1 ,2, • • •  and  then  extrapolates  at  z  *  -®°  .  Schellbach  retains  Stirling's 
approach  casting  the  method  in  an  elegant  algorithmic  form.  The  obvious 
change  of  variable  x  -  rz  transforms  the  problem  into  our  polynomial 
interpolation  in  the  points  of  a  geometric  progression  (Theorem  1  below). 

Quite  recently  I  noticed  the  close  connection  with  the  elegant  algorithm 
of  Romberg  (Theorem  2)  for  which  algorithm  we  refer  to  the  important  paper  [1] 
of  Bauer,  Rutishauser  and  Stiefel.  Also  recently  I  noticed  that  Runge  [3] 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-8Q-C-41 . 


also  solves  the  sane  interpolation  problem  without  stating  this  fact 
explicitly  {Theorem  3).  Rather  he  applies  the  idea  of  the  "Richardson 
deferred  approach  to  the  limit"  20  years  before  Richardson,  working  out  error 
estimates.  We  illustrate  these  connections  with  numerical  examples* 


I.  Stirling  and  Schell bach 


2.  The  Stlrllng-Schellbach  algorithm!  Let 

(2.1)  x^  -  aq*  ,  (k  -  0,1,**#),  a  *  0,  0  <  |q|  <  1 

ba  a  geometric  progression.  Let  f(x)  be  a  function  which  is  analytic  and 
regular  at  the  origin  x  «■  0  /  the  problem  being  to  describe  explicitly 
Newton's  series  (2.1).  Particular  interest  for  applications  is  to  obtain 
approximations  to  the  value  of  f(0)  by  polynomial  extrapolation  at  x  -  0 
from  the  n  +  1  data 

(2.2)  ^  -  f( aqk)  ,  (k  -  0,1, •••  ,n) . 

Following  Schellbach  (5,  §157,  p  280)  we  define  the  q- differences  D*uk 
recursively  by 

(2.3)  0%  “  ^""^k+l  -  ,  (m  -  1,2,“»), 

and  arrange  these  in  a  triangular  array 

u0 

PuQ  -  u,  -  u0 

u1  V\  »  Pu,  -  qPu0 

(2.4)  Pu,  -  Uj  -  u,  P3u0  "  p2“l  “  «2p2u0. 

u2  P2Uj  »  pu2  -  qPuj 

Pu2  -  u3  -  u2 

u3 


Zn  terms  of  these  differences  we  may  state 

1.  For  the  interpolation  points  (2.1)  the  divided  differences  are 


Proof.  Indeed,  for  n  »  1 


S«.*j)  -  -  g,) 

aq  -  a 


Pu„ 


U1  ~  °D _ 

a(q-1)  “  a(q-1 ) 


Assuming  (2*5)  correct  for  n-1  we  have 


-  n-1 , 

f (a,aq,««» ,aq  ) 


nn-1 
V  u. 


n-1  n-t  n-1  .  '  n-1  n-2, 

a  (q  — 1 ) (q  -q) • • • (q  -q  ) 


and  therefore  also 


2  n 

f ( aq, aq  ,«»«,aq  ) 


r»n- 1 
v  u. 


_n-1  n-1,  n-1  4.  ,  n-1  n-2, 

a  q  (q  -1)**»(q  -q  ) 


These  Imply  that 


f(a,aq2,...,aqR)  -  -IteSi.*.*  ligSL-hfi&A * » *  > 

n 

aq  -  a 

«n-1 
V  u. 


1 _  r  1  _  -n-1  i _ 1 

n  n-1  u0*  n-1,  n-1  . , 


,  n  « ,  *  n-1 

a(q  -1 )  q 


.........  .  n-1  n-2 , 

a  (q  -1 )• • • (q  -q  ) 


t>V 


n,  n  . , ,  n  ,  ,  n  n— 1 . 

a  (q  — 1 ) (q  -q)»»»(q  -q  ) 


which  proves  (2.5)  by  induction. 

Using  Newton* 8  expansion  (1.1)  we  immediately  obtain 
Theorem  1.  (Stirling-Schellbach) .  For  the  interpolation  points  (2 
Newton' s  series  (1.1)  becomes 

00 

(2.6)  f(x)  -  l  Pr(x)  PV  , 

n«0 

where 


,4. 


(2.7) 


p  (x)  ,  (x-a)(»-aq)»«»(x-aqn"1_) -  # 

”  a“(qn-1 )(qn-q)***(qn-qn  1 ) 


In  particular,  for  x  -  0  we  obtain  the  expansion 


(2.8)  f(0) 


U0  + 


^0 

i-q 


P2^ 


pV 


(1-q)(1-q  ) 


+•  •  •+ 


n. 


+  •< 


(1_q)...(1_q  ) 


We  propose  to  call  Stlrling-Schellbach  .algorithm  the  sequence  to  sequence 
transformation 

(2.9)  <un>  ♦  (vR) 

which  transforms  the  sequence  (u^)  into  the  sequence  (vn)  of  partial  sums 
of  the  series  (2.8),  hence 


(2.10) 


°U0 

Uq+^-+...  + 


P  u„ 


(1_q)...(1_q  ) 


(n  *  0, !,•••). 


It  is  not  difficult  to  show  that  (2.9)  is  limit  preserving  in  the  sense 

that  if  u  ♦  1  ,  then  also  v  ♦  l  .  However,  the  converse  is  also  true: 
n  n 

If  v  ♦  l  ,  then  also  u  ♦  i.  This  shows  that  (2.9)  can  not  be  used  as  a 
n  n 

limiting  method  that  changes  some  divergent  sequences  (v^)  into  convergent 
ones.  Rather  the  importance  of  (2.9)  lies  in  another  direction:  It  speeds 
up  the  convergence  of  some  slowly  convergent  sequences. 

A  dramatic  example  of  acceleration  of  convergence  was  given  by  Stirling 
himself . 

3.  Stirling's  computation  of  ir  .  We  interpolate  the  entire  function 

2  x/x  1  3 

(3.1)  f(x)  -  —  sin  — —  -  x  -  —  if  x  +  ••• 

/x 

at  the  points 

v  1 

(3.2)  xk  «  q  ,  (k  -  0,1,  — ,n),  where  q-^. 


-5 


Observing  that  the  area  A(m)  of  a  regular  polygon  of  m  sides  inscribed  in 
the  unit  circle  is  given  by 

,  .  .  n  ir 

A(m)  “  m  sin  —  cos  — 

iq  m 


we  find  that 

(3.3)  uk 
and  therefore 

(3.4) 


..,1  .  k+ 1  .  ir _  -k+2  ,  jn _  _rr _ 

f(-)  -  2  sin  —  -  2  sin  —  cos  — 

4  2  2  2 


,  .  1  ,  . ,_k+2. 

-  f(“£>  ”  A(2  >. 

4 


Stirling  computes  the  areas 


A(4 ) ,  A( 8 ) ,  A( 16 ) ,  A(32 ) ,  A(64), 


each  requiring  a  square  root  extraction  to  derive  it  front  the  previous  one. 

We  use  the  Texas  Instruments  SR-50  and  working  with  8  decimal  we  find  that 

(2.4)  becomes 

(3.5) 

uQ  -  A(4 )  -  2  . 

.8284  2712 

u.  -  A(8)  -  2.8284  2712  .0259  3356 

.2330  4034  .0000  9675 

u2  *  A(  16 )  -  3.0614  6746  .00171760  .0000  0006 

.0599  7769  .0000  0157 

u3  -  A(32 )  -  3.1214  4515  .0001  0892 

.0151  0334 

u4  -  A(64 )  -  3.1365  4849 


We  find  the  relevant  q-dif ferences  P^Ug  ,  (k-0,1,2,3,4)  in  the  top 
diagonal  of  this  array.  Forming  the  partial  sums  (2.10)  we  find  the 
approximations 


-6 


ir 


(3.6) 


v1  -  3.1045  6949 


V2  -  3.1414  5277 
V3  -  3.1415  9256 
v4  -  3.1415  9265 

Notice  the  rapid  convergence  of  the  vn  to  ir  ,  the  value  of  v4  having  all 
its  decimals  correct. 

Remarks .  1.  Stirling  [6,  page  133) ,  and  also  Schellbach  who  reproduces 

Stirling's  computations  [5,  page  286}/  uses  also  the  next  area  u5  *  A(128). 
For  some  reason  they  use  in  their  algorithm  these  values  in  their  reverse 
order  A(128),  A(64),  •••,  A(4).  This  reversal  requires  to  replace  q  -  ty4 
by  q  *  4  .  (Extrapolation  at  +*J )  How  does  Schellbach  justify  the  choice 
of  q  «  4?  He  argues  as  follows:  "Since  each  difference  is  nearly  4  times  as 
large  as  the  previous  difference,  the  choice  of  q  »  4  will  result  in  small 
values  of  the  higher  differences  and  lead  to  rapid  convergence".  They  compute 
with  14  decimals  and  obtain  it  with  14  correct  decimals.  Arranging  the 
algorithm  in  the  natural  order  of  (3.5),  and  extending  it  to  include  the  next 
area  u^  •  A(128),  a  fairly  easy  estimate  will  show  the  error  to  be 

| it  -  v5l  <  1.068  x  10’16  . 

Schellbach  devotes  an  entire  chapter  [5,  275-294]  to  Stirling's 
interpolation  series  and  its  applications  and  concludes  the  chapter  by  saying 
"...  this  series,  which  seems  to  have  escaped  so  far  the  attention  of 
mathematicians,  appears  to  be  of  exceptionally  high  practical  importance". 


2.  We  want  to  estimate  the  error  f(0)  -  vn  of  polynomial  interpolation 


of  f( x)  at  x  =  0  .  From  (1.1)  we  know  the  error  at  x  to  be 


f(x 


0' 


,x  ,x)(x-x  )« 
n  u 


(x-XnJ 


(x-xn) •••(x-x  ) 

u  n 


2iri  £ 


f  ( z) dz 


(z-x0). 


(z-x  )(z-x) 
n 


where  T  is  a  closed  curve  containing  x  and  the  interpolation  points.  In 
particular,  for  x  =  0  ,  we  obtain  an  estimate 

n( n+1 ) 

(3.7)  |  f  (0)  -  vj  <  C  •  Kn|q|  2 

valid  for  all  n  .  Here  C  and  K  depend  only  on  f ( x) . 


II.  Romberg 

4.  The  Romberg  algorithm.  We  use  with  slight  mcdi?  : cations  the 
notations  of  the  beautiful  paper  [1]  of  Bauer,  Rutishauser  and  Stiefel.  Let 
r  be  a  constant  such  that 


(4.1) 


\r\  >  1. 


Starting  from  the  column  of  values  we  form  the  Romberg  triangular  array 


(4.2) 

„(0) 


,<0) 


rR^1 '  -  R<0) 
0 _ 0 _ 

r-1 


,(0) 


a^-r^ 

72 .  r 


,(1) 


rR(2)-R(1) 

Ro  Ro 

r-1 


>(0) 


rV1*-^0) 

2  2 

r3  -  1 


(2) 


.(1) 


r2R;2)-R;i) 

/. ; 


.( 2 ) 


rR(3)-R{2) 

rRo  Ro 

r  -  1 


(3) 


the  general  definition  being 


rVk+1>  -  R(k) 

(4.3)  RU)._-1 - -I  or  R(k) 

m  m  ,  m 

r-1 


,(k+1 ) 

m-1 


r-mR(k), 

m-1 


-m 


-9- 


Our  main  result  Is 


Theorem  2.  The  Stlrling-Schellbach  algorithm  (2.10)  is  equivalent  with 
the  Romberg  algorithm  (4.2)  such  that 

(4.4) 

This  means  the  following:  If  we  identify  the  first  columns  of  the  arrays 

(2.4)  and  (4.2)  ,  so  that 

(4.5) 

then  the  elements  of  the  leading  diagonal  of  (4.2)  are  respectively  equal  to 


1 

r  =  — 

q 


u  =  R*m>  ,  (m  =  0,1,«»»)» 
m  u 


the  partial  sums  of  the  series  (2.8),  l.e. 

(4.6)  v  =  R(0>  ,  (m  =  0, 1, •••). 

m  m 

We  give  two  proofs. 

First  proof.  My  colleague  C.  de  Boor  pointed  out  to  me  the  following 
remark :  If  we  apply  the  Neville  algorithm  for  a  geometric  progression  xm  « 
aq”1  and  for  interpolation  at  x  =  0  ,  then  Neville's  fractions  simplify  and 
become  identical  with  the  elements  of  Romberg's  algorithm.  This  connection 
was  already  mentioned  in  (4,  301-302]. 

Second  proof.  This  proof  is  more  direct  but  longer.  The  equation  (4.6) 
is  true  by  definition  if  m  =  0  ,  because  vQ  *  Uq  *  Rq^0^.  To  establish  the 
equation 


Vuo 

u„  +  - —  +  »••  + 

0  1-q 


Pm 


.(0) 


(1_q)...(1_qm)  m 


for  all  m  ,  we  assume  that  it  holds  for  m  -  1  ,  and  we  are  to  show  that 


Vmu 


(4.7) 


( 1 -q ) • • • (1-qm) 


r<0)  -  r(o;  ,  with  r<;>  -  o 

m  m-1  -1 


Let  us  prove  this  by  induction.  For  m  =  0  this  reduces  to  uQ  *  R^ 


(0) 


Assuming  (4.7)  true  for  m-1  we  have 


(4.8) 


P*"1 


( 1-q) • • • ( 1-q  ) 


-  R(0>  -  P(0) 
nr-V  Knr1  m-2  ' 


and  we  are  to  derive  from  it  that  also  (4.7)  holds,  i.e.. 


(4.9) 


flr1  n 

V  u1  -  q 


■  1  m-1 

V 


(1-q) •••(1-q  ) 


,(0)  „(0) 


n 


-  R 


m-1 


However,  the  assumption  (4.8)  also  implies  that 


,m— 1 


(1-q)... (1-q  ) 


R110  -  R(lt) 
1.  m-1  m-2 


and  in  particular,  for  k  «*  1,  that 


n  m-1 
V  u. 


(4.10) 


( 1-q) • • • ( 1-q  ) 


R(1)  -  R(1) 
m-1  m-2 


Using  (4.8)  and  (4.10),  (4.9)  becomes 


(4.11) 


R(1!  -  R(1i 
m-1  m-2 


1  -  q 


m 


-  q 


m-1 


R(0)  -  R<0) 

m-1  m-2 

1  -  q 


R<°>  -  R(0> 
m  m-1 


However,  the  second  equation  (4.3),  with  r”1  *  q  ,  shows  that 


(1)  -  (1  -  qm)R(0)  -  -n»(0) 


m-1 


*  "  +  q  Rm-1 


c(1)  M  „®-1,B(°)  *  ►1B(°) 

Rm-2“n_q  Rm-1  +  q  Rm-2 


Substituting  into  (4.11)  we  obtain 


11 


(1-q  )Rn  +  q^a,1  -  (1-q  -  q 

1  -  q 


-  q 


,  *<0] 

■1  m-1 


-  R 


(0) 


m-2 


1  -  q 


m 


r(o)  -  r(o; 

m  m-1 


which  sinplifies  to 


"-X01*  X°l  - *, 


(0) 

m-1 


1  -  q" 


R(o)  -  r(o; 
m  m-1 


which  is  visibly  an  identity • 

Remarks .  1*  Even  though  the  two  algorithms  (2.10)  and  (4.2) ,  with 

u  “  and  r  -  q-1  ,  solve  the  same  interpolation  problem,  it  is  clear 

that  the  elegant  Romberg  algorithm  is  much  to  be  preferred.  We  illustrate 
this  by  returning  to 

2.  The  computation  of  r  .  With  r  »  q“*  *4  and  for  the  uk  of 
(3.5)  Romberg1 s  triangular  array  becomes 


3.1415  8294 

u4  *  3.1365  4849  . 


We  recognise  in  the  top-diagonal  of  (4.12)  the  values  (3.6),  except  for  some 


rounding  errors,  as  guaranteed  by  Theorem  2 


ZZZ.  Rung* 

5.  Runge1 s  first  problem.  Without  knowledge  of  the  work  of  Stirling  end 
Schellbach,  Runge  considers  in  (3]  the  following  problems.  Let 

(5.1)  f(x)  *  aQ  +  a,|X  +  •••  +  a^x*  +  *  ( |x|  <  k) 

be  regular  in  the  circle  |x|  <  r  .  Let  q  and  a  be  constants  such  that 

(5.2)  0  <  |q|  <  1  ,  0  <  | a |  <  r  . 

Runge* s  Problem  1  is  to  approximate  to  ig  ■  f(0)  in  terms  of  the  n+1 
values 


(5.3)  u^  «  f(aqK)  ,  (k  -  0,1, ••*,n)  . 

This  problem  was  solved  by  Theorem  1  by  the  Stirling-Schellbach  algorithm 
(2.10),  and  also  by  Theorem  2  by  the  Romberg  algorithm  with 

(5.4)  r  -  -  . 

q 

However,  nowhere  does  Runge  mention  this  polynomial  extrapolation  approach. 
Rather  he  proceeds  directly  as  follows,  considering  only  the  case  when 
q  m  V2  •  We  write  the  (5.3)  as 


n-k.  _  nk.  2  2n  2k  , 

f(aq  )  *  aQ  +  a^q  r  +  a2a  q  r  + 


(5.5) 

J  n  n2  nk  n+1  (n+1)n  (n+1)k  .  ..  A  .  . 

+  *na  q  r  "-Vi*  q  r  +  '  (k  -  0,1,...,n). 

He  multiplies  this  equation  with  and  sums  over  all  k  ,  the  objective 

being  to  so  choose  the  as  to  anihilate  the  n  terms  in  a®  (s  - 

1,2,  ***,n).  This  is  seen  to  be  achieved,  provided  the  are  the 

coefficients  of  the  polynomial 

(5.6) 

defined  by  the  conditions 

(5.7) 

Zn  this  way  Runge  obtains  from  (5.5)  the  equation 


^  (x,  .  c  +  c  x  +  •  ••  +  c  x"  -  (r-x)(r2-x)»y(rn-x) 

n'  0  1  n  ,  . , ,  2  4 .  ,  n  . . 

(r-1)(r  -1)»»»(r  -1) 


w  (r)  -  w  ( r2)  -  •••  -  *  (rn)  -  0  ,  C 1 )  -  1  . 
n  n  n  n 


n 


(5.8)  -  l  Cvf(aq  K)  -  afl  +  a^a  q 


n  “  k 


I  g  IIT  |  . 

V»(r  )  +  ••• 
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.  .  n+1 ,  .  .  .n  1+2+***+n  ,  ,  .n  -n(n+1)/2  __  ..  .  _  _ , 

Sine*  <fi(x  )  »  (-1 )  r  ■  (-1)  q  w*  that  Range's 

n 

approximation  *n  to  *0  gives  tha  same  order  of  approximation  (3.7)  a a  tha 


Stirling-Schellbach  vR  .  Tha  following  result  should  therefore  coma  as  no 
great  surprise. 

Thereom  3.  Rungs* s  approximation  &n  of  P(0)  -  Sq  ,  defined  by  (5.8), 
(5.6),  (5.7),  is  identical  with  the  result  of  polynomial  extrapolation,  hence 

(5.9)  «n  *  »»- 

A  proof  of  this  is  fairly  easy  if  we  define  the  interpolated  value  v„ 
by  Lagrange's  interpolation  formula,  and  also  use  the  Gaussian  identity 


(x-r)(x-r2)*«*(x-rn)  -  l  (-1)V[  "  ]r 

0 


V(V+1) 


n-v 


where 


,  n  ...  n-l  ..  .  n-v+i 

r  n  i  (r  — 1 ) (r  -1)«««(r  -1 ) 

l  u  J  v  V—  1 

(r  -1 )(r  -i)...(r-l) 


We  omit  the  details. 

Remarks.  Besides  the  Stirling  computation  of  «  pf  §3  ,  further 
attractive  examples  are  provided  by  the  following  functions. 

1.  The  function  1 

(5.10)  f(x)  -  (1  +  x)x  •  *  +  a^x  +  •••  . 

To  determine  f(0)  ■  •  w*  can  us*  Romberg's  algorithm  choosing  e.g. 
a  ■  1/8,  q  ■  1/16,  hence  r  ■  16.  Notice  that  xk  -  aq*  ■  2_3~4k  and 
therefor*  the  computation  of  uk  ■  f(aqk)  from  (5.10),  requires  only 
successive  squaring. 

2.  The  entire  function 

(5.11)  f(x)  -  “<2X-1)  -  log2  +  a^  ♦  ••• 

will  give  approximations  to  log  2.  with  a  -  1  and  q  •  ^  ,  the  computation 


14< 


of  uk  -  f(aqk)  •  f(1/2k)  requires  only  successive  square  root  axtractions. 
Zn  (5.11)  we  nay  roplaca  2  by  any  integer  H  . 


3.  The  case  when  f(x)  is  an  even  function  of  x  .  If 

(5.12)  f(x)  “  Xq  +  ajX2  +  a4x4  ♦  ••• 

we  define 

2 

(5.13)  g(x)  -  f(/x)  -  aQ  +  a2x  +  a2x  +  ••• 
and  observe  that 

(5.14)  g(a2q2k)  -  f(aqk)  ,  (k  -  0,1,»««). 

Since  g(x)  is  also  regular  at  x  ■  0  ,  we  obtain 

Theoree  4.  If  the  function  (5.1)  is  even*  then  we  nay  apply  to  the 
data  (5.3)  the  aneberq  algorithm  with  q  replaced  by  q2  •  This  eeans  that 
in  Rosibora* s  aloorlthst  (4.2)  we  keep  the  first  oolasne 

-  u,  -  f(aq") 

fixed  and  pass  fro*  r  -  q"1  to  r2  -  q"2  . 

Replacing  q  by  q2  will  clearly  accelerate  the  convergence  of  the 
vn  to  f(0)  -  .  This  important  device  was  already  used  in  §3i  Instead  of 

interpolating  the  even  function 

•  sin  ^  at  x  -  1,  -j  »  -tj  #  ••• 

we  interpolated  the  entire  function  (3.1)  at  x  ■  1,  1/4.  1/42  ,  •••  •  In 

{6  we  shall  again  use  this  device  to  good  advantage. 
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Let 


'a  second  problem:  Computing  the  invars*  functio 


(6.1) 


-  «x)  ■  «lq  +  a1x  +  ajX2  +  ••• 


(  |x |  <  r) 


be  such  that 

(6.2)  w*  can  compute  the  value  of  f(xq)  in  terms  of  f(x)  , 
and  therefore,  successively, 

(6.3) we  can  compute  f(xq*)  (k  ■  0, !,*••)  in  terms  of  y  »  f(x). 
Throughout  this  section  we  assume  that 


(6.4) 


q  ■!/$* 


The  formulae 


,  x 
sin  - 


1  -  cos  x 
2 


r1  -  r  1  -  sin  x 
2 


show  that  sin  x  satisfies  (6.3)  with  q  *  ^  . 

Problem  2.  To  compute  x  _if  y  -  f(x)  is  prescribed. 

Solution.  Rungs  reduces  this  problem  to  Problem  1  of  $5  as  follows. 


Writing 


f ( xt )  -  *q  +  a^xt  +  a2x2t2  +  a3x3t3  +  •••  ,  (a^  *  0)  , 


we  define  the  new  function 


(6.5) 


,|Xt|-*0  *2  2.  *3  3.2 

— - - ■  x  +  —  xt  +  —  Xt  +  ••• 

V  *1  *1 


which  depends  also  on  x  ,  and  for  which 


.  k.  f(wl  }  "  *0  *2  2  k  *3  3  2k  . 

-  9(q  )  •  - £ -  -  x  ♦  —  xq  +  pxq  +•••,  (k  -  0 , 1 ,  *  •  • ) , 

•,q  1  1 


»y  our  assusiptlon  (6.3)  all  these  values  can  be  computed,  and  (6.5)  shows  that 


(6.7) 


*  -  g(0) 
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can  be  obtained  as  the  solution  of  Problem  1  for  the  data  . 

Following  Runge  [3/  222-223]  and  also  Schellbach  [5,  88-90]  we 
illustrate  this  procedure  by 

The  computation  of  an  incomplete  elliptic  integral  of  the  first  kind. 

Let 

(6.8)  x  -  JJ  —  ■  -  -  -  ,  (0  <  c  <  1  ,  0  <  y  <  1/2)  , 


x  -  /y 
0 


/l  -  c2sin26 


where  y  is  prescribed  and  we  are  to  compute  x  .  A  result  of  Legendre  [2, 
vol.  1,  521#  25-26]  is  as  follows:  If  we  determine  acute  angles  t  and  y^ 
from  the  equations 

(6.9)  siny  ■  c  sin  y#  sin  y1  “  sin  ^  /  cos  ^  , 

then 

d6 


(6.10) 


f  ■  /y' 


0  A  -  C2.ln2e 


This  shows  how  the  value  of  the  integral  can  be  halved.  This  operation  can  be 
repeated:  He  determine  successively  angles  Yr_1  and  yn  from  the  equations 


•1»  v,  •  c  Vf 


1  i  /  Yn-1 

sin  yn  -  sin  ~ —  /cos  ^ — 


(6.11) 


(n  -  1 ,2, • •• ,k  t 


y0  -  y0"y)/ 


to  obtain  yfc  such  that 


(6.12) 


_x_  f\ 

jt  "  * 

2  0 


d6 


/l  -  c2sin20 


He  now  invert  the  relationship  (6.8)  to  obtain 
(6.13)  y  -  f(x)  »  x  +  a3x3  +  a5x5  ♦  •••  , 

observing  that  it  is  an  odd  function  which  is  usually  denoted  by  f(x)  *  sn  x. 


Clearly  we  may  rewrite  (6.12)  as 


•  ->lr 


(6.14) 


In  terms  of  (6.13)  our  function  (6.5)  becomes 


(6.15) 


„  x  +  a  +  a  x^t4  +  ••• 

t  3  5 


and  in  particular  (6.6)  may  be  written  as 


,  1  .  _k. ,  x  t  .  k 

uk  ’  “ 2  f(pr>  ’ 2  yk 


(6.16) 


x  +  a3x32~2k  +  a5x52_4k  +  •••  ,  (k  -  0, !,•••)  . 


This  equation  shows  that  u^*  x  as  k  ♦  •  ,  and  u^  is  Legendre's 
approximation  to  x  .  However,  this  approximation  can  be  much  improved  by  our 
extrapolation  procedure. 

A  numerical  example.  Let  us  evaluate  the  incomplete  elliptic  integral 


(6.17)  x  -  J4  . .  .  where  c  -  !^  and  y  -  j  . 

0  /l  -  V4  sin20 

From  Legendre's  Table  VIII  [2,  vol.  3,  page  339]  we  obtain  the  12  place  value 


(6.18) 


x  *  .80436  61012  32. 


The  following  computations  were  made  in  double  precision  by  Fred  W.  Sauer,  of 
the  MRC  Computing  Staff. 

We  solve  the  equations  (6.11)  with  k  -  4  and  find  the  following  angles 
in  radians 


Y  -  .36136  71239  06708 


Y 1  -  .19575  59443  11023 

Y_  -  .09987  08442  57596 
2 

Y3  -  .05018  82616  25625 


y1  -  .39956  34259  31126 

y„  -  .20075  55982  73493 

y3  -  .10050  35008  99138 

yA  -  .05026  75900  94179  . 
4 


From  (6.16)  we  get  the  values 


y  -  .78539  81633  97448 


u0  * 

u1  -  2y1  -  .79912  68518  62251 

(6.19)  u2  -  4y2  -  .80302  23930  93970 

u3  -  8y3  -  .80402  80071  93103 

U4  •  16y4  -  .80428  14415  06865. 

The  last  value  u4  Is  Legengre's  3-place  approximation  to  x  .  We  can  now 
improve  this  approximation  in  two  different  ways* 

1.  The  Romberg  algorithm  applied  to  the  data  (6.19)  with 

(6.20)  r  -  2  gives  the  approximation  R4^  “  .80436  56250 
which  is  seen  to  be  correct  to  6  decimal  places. 

2.  Because  the  g(t)  of  (6.15)  is  an  even  function  we  may  use  Theorem  4 
and  find  that  the  Romberg  algorithm  with 

(6.2*1 )  r  -  22  -  4  gives  the  approximation  R^°*-  .80436  61012  29163 

which  is  seen  to  be  correct  to  nearly  12  decimals. 

To  conclude  it  seems  worthwhile  to  recall  that  Romberg  invented  his 
algorithm,  also  with  r  ■  4  ,  for  the  evaluation  of  a  definite  integral 
(6.17)  in  terms  of  its  binary  trapezoidal  sums  (See  [1]).  Writing 
F(0)  *  (1  -  V4  sin28)  ^2  ,  we  consider  the  sums 

2^- 1 

T  -  J  2~k{  ^F(O)  +  l  «  1  2'k)  +  V(J)1  -  (k«0,1,*«M. 

2*  8-1 

Dividing  [0,j]  into  24  -  16  equal  parts,  we  can  compute  the  5  sums 

(6.22)  T^,  T2,  T4,  Tg,  T^g  . 

Zf  we  now  apply  to  the  column  of  the  sums  (6.22)  the  Romberg  algorithm  with 
r  -  4  we  find  with  double  precision  the  approximation 

(6.23)  R<0)  -  .80436  61012  31069. 

4 

A  comparison  with  Legendre's  value  (6.18)  shows  that  (6.23)  is  even  slightly 
more  accurate  than  (6.21). 
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